%% Generate figures
clear;

% Settings
[bParams, oParams, setts]=model_settings(2);
dir = setts.dir;
pbar = oParams.pbar;
fsep = filesep;
products = (0:pbar)';
bin_n = 10;
bin = (1:bin_n)';

% Load results
load(os_dir(setts.dir,'\cs_table.mat'));
load(os_dir(setts.dir,'\nodebt_results.mat'));
load(os_dir(setts.dir,'\static_results.mat'));

setts.dir = dir;

%% Figure 7 (firm size plots)

[N_nd,edges] = hist(dataSim_nd.pt(:),1:pbar);
distribution_nd = [0 N_nd/length(dataSim_nd.pt(:))]';
[N_ds,edges] = hist(dataSim_ds.pt(:),1:pbar);
distribution_ds = [0 N_ds/length(dataSim_ds.pt(:))]';

table_size = table(products,distribution_nd,distribution_ds);
writetable(table_size,os_dir(setts.dir,'data\figure_7.txt'),'Delimiter',' ');

%% Figure 6 (investment plots)

[bParams, oParams, setts] = model_settings(2);
p_lambda_ds_fb = [NaN; max(underinvestment(oParams,bParams,setts,f_ds,dataSim_ds.pt,dataSim_ds.ct,dataSim_ds.lambda),0)'];
[bParams, oParams, setts] = model_settings(2);
bParams.beta = 4.503;
p_lambda_ds_fb_betai_high = [NaN; max(underinvestment(oParams,bParams,setts,f_betai_high,dataSim_betai_high.pt,dataSim_betai_high.ct,dataSim_betai_high.lambda),0)'];
[bParams, oParams, setts] = model_settings(2);
bParams.H = 4.138;
p_lambda_ds_fb_H_high = [NaN; max(underinvestment(oParams,bParams,setts,f_H_high,dataSim_H_high.pt,dataSim_H_high.ct,dataSim_H_high.lambda),0)'];

table_fb = table(products,p_lambda_ds_fb,p_lambda_ds_fb_betai_high,p_lambda_ds_fb_H_high);
writetable(table_fb,os_dir(setts.dir,'data\figure_6.txt'),'Delimiter',' ');

%% Figure 8 and 9 (f_entrant plots)

list = {'beta', 'gamma', 'pi', 'H'};

for i=1:length(list)
	eval([list{i},'=readtable(os_dir(setts.dir,''\data\debt\',list{i},'.txt''));']);
	eval([list{i},'_nd=readtable(os_dir(setts.dir,''data\nodebt\',list{i},'.txt''));']);
	eval([list{i},'_vec=table2array(',list{i},'(:,1));']);
	eval(['f_entrant=table2array(',list{i},'(:,3));']);
	eval(['f_entrant_nd=table2array(',list{i},'_nd(:,3));']);
	eval(['f=table2array(',list{i},'(:,2));']);
	eval(['f_nd=table2array(',list{i},'_nd(:,2));']);
	eval(['m_entrant=table2array(',list{i},'(:,10));']);
	eval(['m_entrant_nd=table2array(',list{i},'_nd(:,10));']);
	eval(['m_incumbent=table2array(',list{i},'(:,9));']);
	eval(['m_incumbent_nd=table2array(',list{i},'_nd(:,9));']);
	f_incumbent=f-f_entrant;
	f_incumbent_nd=f_nd-f_entrant_nd;
	products_incumbent=(f_incumbent)./m_incumbent;
	products_incumbent_nd=(f_incumbent_nd)./m_incumbent_nd;
	products_entrant=(f_entrant)./m_entrant;
	products_entrant_nd=(f_entrant_nd)./m_entrant_nd;
	change_products_incumbent=(products_incumbent-products_incumbent_nd)./products_incumbent_nd;
	change_mass_incumbent=(m_incumbent-m_incumbent_nd)./m_incumbent_nd;
	change_incumbent=(f_incumbent-f_incumbent_nd)./f_incumbent_nd;
	change_products_entrant=(products_entrant-products_entrant_nd)./products_entrant_nd;
	change_mass_entrant=(m_entrant-m_entrant_nd)./m_entrant_nd;
	change_entrant=(f_entrant-f_entrant_nd)./f_entrant_nd;
	share_entrant=(f_entrant-f_entrant_nd)./(f-f_nd);

	diff_f_incumbent=f_incumbent-f_incumbent_nd;
	diff_f_entrant=f_entrant-f_entrant_nd;
	diff_m_incumbent=m_incumbent-m_incumbent_nd;
	diff_m_entrant=m_entrant-m_entrant_nd;
	diff_products_incumbent=products_incumbent-products_incumbent_nd;
	diff_products_entrant=products_entrant-products_entrant_nd;

	eval(['table_fent = table(',list{i},'_vec,change_entrant,share_entrant,change_mass_entrant,change_products_entrant,change_products_incumbent,change_mass_incumbent,change_incumbent,diff_f_incumbent,diff_f_entrant,diff_m_incumbent,diff_m_entrant,diff_products_incumbent,diff_products_entrant);']);
	eval(['writetable(table_fent,os_dir(setts.dir,''data\figure_8_and_9_',list{i},'.txt''),''Delimiter'','' '');']);
end

%% Figure 5 (entrant value vs f)

lowerb = Mom_ds{:,1}*0.25;
upperb = Mom_ds{:,1}*2.5;
n_points = 100;

% f grid, static debt
[bParams, oParams, setts] = model_settings(2);
f_grid_ds = linspace(lowerb,upperb,n_points)';
upMat = updatematrix(bParams,oParams,setts);
entrant_vs_f_ds = zeros(size(f_grid_ds));

for i=1:length(f_grid_ds)
    entrant_vs_f_ds(i) = entrant_value(bParams,oParams,setts,upMat,f_grid_ds(i));
end

% f grid, static debt, low beta
[bParams, oParams, setts] = model_settings(2);
bParams.beta = 3.2528;
f_grid_ds_beta_low = linspace(lowerb,upperb,n_points)';
upMat = updatematrix(bParams,oParams,setts);
entrant_vs_f_ds_beta_low = zeros(size(f_grid_ds_beta_low));

for i=1:length(f_grid_ds_beta_low)
    entrant_vs_f_ds_beta_low(i) = entrant_value(bParams,oParams,setts,upMat,f_grid_ds_beta_low(i));
end

% f grid, static debt, n high
[bParams, oParams, setts] = model_settings(2);
oParams.n = 3;
f_grid_ds_n_high = linspace(lowerb,upperb,n_points)';
upMat = updatematrix(bParams,oParams,setts);
entrant_vs_f_ds_n_high = zeros(size(f_grid_ds_n_high));

for i=1:length(f_grid_ds_n_high)
    entrant_vs_f_ds_n_high(i) = entrant_value(bParams,oParams,setts,upMat,f_grid_ds_n_high(i));
end

table_f_entrant = table(f_grid_ds, entrant_vs_f_ds, f_grid_ds_beta_low, entrant_vs_f_ds_beta_low, f_grid_ds_n_high, entrant_vs_f_ds_n_high);
writetable(table_f_entrant,os_dir(setts.dir,'data\figure_5.txt'),'Delimiter',' ');

%% Figure 3 (PD graphs)

p0=[0:15]';
p0_max=15;
pds_ds = pd_plot(equity_mat_ds,optimal_coupon_ds,p0_max,setts);
[bParams, oParams, setts] = model_settings(2);
bParams.beta = bParams.beta*0.8525;
[Mom_betai_pd, dataSim_betai_low_low, f_betai_pd, ~, ~, equity_mat_betai_pd, ~, ~, ~, ~, optimal_coupon_betai_pd] = moments(bParams,oParams,setts);
pds_ds_betai = pd_plot(equity_mat_betai_pd,optimal_coupon_betai_pd,p0_max,setts);
[bParams, oParams, setts] = model_settings(2);
oParams.pi = 0.25;
[Mom_pi_pd, dataSim_pi_pd, f_pi_pd, ~, ~, equity_mat_pi_pd, ~, ~, ~, ~, optimal_coupon_pi_pd] = moments(bParams,oParams,setts);
pds_ds_pi = pd_plot(equity_mat_pi_pd,optimal_coupon_pi_pd,p0_max,setts);
table_pd = table(p0,pds_ds,pds_ds_betai,pds_ds_pi);
writetable(table_pd,os_dir(setts.dir,'data\figure_3.txt'),'Delimiter',' ');

%% Figure 10 (elasticity figure)

H_ds=readtable(os_dir(setts.dir,'data\debt\H_ext.txt'));
H_vec_ds=table2array(H_ds(:,1));
exit_rate_H_ds=table2array(H_ds(:,8));
M_i_H_ds=table2array(H_ds(:,9));
f_H_ds=table2array(H_ds(:,2));
M_e_H_ds=table2array(H_ds(:,10));
f_e_H_ds=table2array(H_ds(:,3));

[bParams, oParams, setts] = model_settings(1);
[Mom_nd, dataSim_nd, f_nd, entrant_nd, lambda_e_nd, equity_mat_nd, debt_mat_nd, lambda_mat_nd, coupon_mat_nd, optimal_firm_nd, optimal_coupon_nd] = moments(bParams,oParams,setts);
M_i_nd = table2array(Mom_nd(:,8));
M_e_nd = table2array(Mom_nd(:,9));
exit_rate_nd = table2array(Mom_nd(:,7));
H_nd = bParams.H;

[bParams, oParams, setts] = model_settings(2);
elasticity_ds_scaled = log((M_i_H_ds(2:end) .* exit_rate_H_ds(2:end))/(M_i_nd .* exit_rate_nd))./log(H_vec_ds(2:end)/H_nd);
plot(elasticity_ds_scaled,f_H_ds(1:end-1),elasticity_ds_scaled,f_nd*ones(length(elasticity_ds_scaled),1))
elasticity_plot = elasticity_ds_scaled(1:end);
f_ds_plot = f_H_ds(1:end-1);
f_i_ds_plot = table2array(H_ds(1:end-1,4));
f_nd_plot = f_nd*ones(length(elasticity_ds_scaled),1);
f_i_nd_plot = table2array(Mom_nd(:,3))*ones(length(f_i_ds_plot),1);
f_e_nd_plot = table2array(Mom_nd(:,2))*ones(length(elasticity_ds_scaled),1);
f_e_ds_plot = table2array(H_ds(1:end-1,3));
m_i_ds_plot = M_i_H_ds(1:end-1);
m_i_nd_plot = M_i_nd*ones(length(elasticity_ds_scaled),1);
m_e_ds_plot = M_e_H_ds(1:end-1);
m_e_nd_plot = M_e_nd*ones(length(elasticity_ds_scaled),1);
prod_i_ds_plot = f_i_ds_plot./m_i_ds_plot;
prod_i_nd_plot = f_i_nd_plot./m_i_nd_plot;
prod_e_ds_plot = f_e_ds_plot./m_e_ds_plot;
prod_e_nd_plot = f_e_nd_plot./m_e_nd_plot;

table_p = table(elasticity_plot, f_ds_plot, f_i_ds_plot, f_e_ds_plot, f_nd_plot, f_i_nd_plot, f_e_nd_plot, m_i_ds_plot, m_i_nd_plot, m_e_ds_plot, m_e_nd_plot, prod_i_ds_plot, prod_i_nd_plot, prod_e_ds_plot, prod_e_nd_plot);

change_entrant = (table_p.f_e_ds_plot - table_p.f_e_nd_plot)./table_p.f_e_nd_plot;
change_incumbent = (table_p.f_i_ds_plot - table_p.f_i_nd_plot)./table_p.f_i_nd_plot;
change_mass_entrant = (table_p.m_e_ds_plot - table_p.m_e_nd_plot)./table_p.m_e_nd_plot;
change_mass_incumbent = (table_p.m_i_ds_plot - table_p.m_i_nd_plot)./table_p.m_i_nd_plot;
change_products_entrant = (table_p.prod_e_ds_plot - table_p.prod_e_nd_plot)./table_p.prod_e_nd_plot;
change_products_incumbent = (table_p.prod_i_ds_plot - table_p.prod_i_nd_plot)./table_p.prod_i_nd_plot;
elasticity = table_p.elasticity_plot;

table_el = table(elasticity, change_entrant, change_incumbent, change_mass_entrant, change_mass_incumbent, change_products_entrant, change_products_incumbent);
writetable(table_el,os_dir(setts.dir,'data\figure_10.txt'),'Delimiter',' ');

%% Smoothing figures

run(os_dir(pwd,'\functions results\smoothing_cs.m'))